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This thesis is concerned with a comparative study of discrete 
time filters using the theories of Wiener-Kolmogorov , Bode-Shannon, 
and Kalman, applied to the filtering of a non-stationary random 
signal in the presence of measurement noise. Programs are developed 
for the simulation of these systems and signals on a digital computer. 
Their filtering properties are compared for a random input signal 
with known steady state characteristics, starting at initial time 
t = t . The results show that when the gain of the Kalman filter 
is equal to its steady state value, the Kalman and Wiener-Kolmogorov 
filters perform identically. For large initial errors the Kalman 
filter, with large initial gain, gives the best transient response; 
for small initial errors the Kalman and Wiener-Kolmogorov filters 
are essentially equivalent in their transient responses. 
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CHAPTER I 



INTRODUCTION 
1. OBJECTIVE 

This thesis examines and compares several techniques for the 
filtering of random signals in the presence of measurement noise. 

In diagram form, the problem is as shown in Figure 1; x(t) is the 
random signal to be estimated, v(t) is a random measurement noise, 
£(t) is the estimated or filtered signal, and e(t) is the error in 
estimation. The criteria for the optimum filter is one which 
reduces the mean square error to a minimum. The following techniques 
are studied and compared with a common example: 

(1) Discrete Wiener-Kolmogorov filter 

(2) Discrete Kalman filter 

(3) Bode-Shannon discrete time filter 

For all of the filters it is assumed that the statistical 
properties of the signal and the noise are known by their correlation 
functions, or their power density spectra, or by equivalent white 
noise excited dynamic models. It is also assumed that the signal 
and the measurement noise are uncorrelated. Wiener-Kolmogorov theory 
assumes that the signal and the noise are stationary; that is, their 
statistical properties do not vary with time. This implies that all 
signals are initiated at t = Kalman and Bode-Shannon filtering, 

on the other hand, assume that signals start at some initial time, 
t^, and that as time progresses the stationary probability charac- 
teristics of the signals are established. In addition, Kalman (as 



9 




10 



BASIC FILTER PROBLEM 



well as Bode-Shannon) assumes that the random signal process is 

Markovian, that is, can be thought of as generated by a linear 

* 

dynamic system excited by white noise. During the transient 
stage, from t = t^ until the statistical properties are established, 
the Kalman and Bode-Shannon filters behave as time-varying devices. 
In the steady state, all of the filters can be expected to behave 
identically since they each establish a minimum mean square error 
estimate . 

In order to test and compare the operation of these filters, 
particularly for a non- stationary input, the Wiener-Kolmogorov and 
Kalman filters have been designed for discrete time simulation on 
a digital computer for filtering the same random input signal and 
measurement noise. The output of each filter and the actual random 
input signal have been compared starting at time t = t^ until steady 
state conditions are established. 

The type of filtering problem studied can be compared with a 
basic tracking problem where a target, with statiscally known 
maneuvering properties, suddenly appears. It is desired to track 
the target and, in particular, to investigate how long it takes the 
output to reach the steady state minimum mean squared error. The 
results of this investigation are presented in this thesis together 
with a summary of pertinent background theory. In addition, the 
discrete time formulations for the Wiener-Kolmogorov theory have 
been developed and programmed for the specific example treated. 

*R. E. Kalman, New Methods and Results in Linear Prediction 
and Estimation Theory , Technical Report No. 61-1, (Baltimore, 
Maryland: Research Institute for Advanced Study, 1961) p. 3. 
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The basic theory for discrete Kalman and Bode-Shannon filtering 
appears in the literature and has been adapted to the example. 
Programs necessary for generating the stochastic signals from 
white noise have also been developed. 

2. HISTORICAL BACKGROUND 

The many present day theories for the filtering of time 
sequences are basically derived from the original work of Wiener 
and Kolmogorov published in 1941 and 1942 ? 9 ^ 5 ^ They obtained 
the direct solution of the filter problem independently and 
simultaneously by the use of the calculus of variations. This 
approach led to the filter's impulse response being characterized 
by the solution of a difficult integral equation known as the 
Wiener-Hopf equation. 

One of the next major developments in the solution of the 

2 17 

filtering problem was accomplished by Bode and Shannon 5 in 1950 
with the simplification of the derivation of the Wiener-Kolmogorov 
filter. They used frequency domain analysis and conventional 
circuit theory concepts to interpret mathematical operations in 
physically intuitive terms. Bode and Shannon also applied this 
frequency domain technique to discrete signals with constant 
sampling intervals . 

In most cases, the solution of the Wiener-Hopf integral equation 

is a formidable task. With the advent of the digital computer, an 

alternate approach was sought to avoid this difficulty. In 1960 
9 10 11 11 

Kalman 5 ’ and Bucy attacked the filtering problem from the 

state variable point of view. They evolved a set of differential 
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equations, equivalent to the Wiener-Hopf integral equation, whose 
solution yields the optimum filter. These differential equations 
may be synthesized in a sequential fashion and therefore digital 
computer solutions easily obtained. 



13 



CHAPTER II 



THE WIENER- KOLMOGOROV FILTER 



There are several methods for deriving the continuous Wiener- 
Kolmogorov filter. Appendix A summarizes the solution which leads 
to the Wiener-Hopf integral equation.^ 



cp (t) 

T zx 






(T) CD (t - T)dT, t S 0 
zz 



( 2 . 1 ) 



H(t) is the desired filter impulse response, cp (t) is' the auto- 

z z 

correlation function of the signal plus noise, and cp (t) is the 

z x 

crosscorrelation function of the signal plus noise, and the signal. 

The requirements that must be fulfilled to use Eq. 2.1 are 
(1) that all signals are stationary and exist for all time; (2) that 
the signal and the noise are uncorrelated. In order to solve the 
equation, the autocorrelation function of the signal and the auto- 
correlation function of the noise must be known. These may be 
obtained by taking the Fourier transform of the respective power 

■>v 

spectral densities, as given by the Wiener-Khintchine equations: 



W (to) = 
x 



h / Cf, xx (T)e ' : ' l " TdT - 






cp (t) 
T xx v 7 



George C. Newton, Jr., Leonard A. Gould, and James F. Kaiser, 
Analytical Design of Linear Feedback Controls (New York: John Wiley 

and Sons, Inc., 1957) p. 144. 

i< 

John L. Stewart, Fundamentals of Signal Theory (New York: 
McGraw-Hill Book Company, Inc., 1960) p. 306. 
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cp (t) = 
T xx 



f W (U))e + ^ T d<JU = 2 ttj« 1 


W (U>) 


L x 


S' 



where W (o>) is the power density spectrum of x(t). Even when 

x 

sufficient knowledge of the signal is available, the solution of 
Eq. 2.1 is generally an extremely difficult task. 

1. SOLUTION BY SPECTRUM FACTORIZATION 

The solution of Eq. 2.1 using the technique known as spectrum 

2 1 17 

factorization has been developed by several authors. 5 ’ The 

autocorrelation function of the signal plus noise, cp (s), is 

z z 

factored so that 



* M (s) ■ V zz (s) 2 z (s) 



( 2 . 1 - 1 ) 



where cp (s) contains those factors of cp (s) which have poles and 
z z z z 

zeros in the left half plane. The term 



cp (s) 
T zx 

T zz 



(2.1- 2a) 



is defined by those terms of the partial fraction expansion of: 



cp (s) 
T zx 






which have poles in the left half of the s-plane. That is 



9 



zx 



9 



zz 



cp (s) 




ZX 


+ 


cp (s) 

T zz 






+ 






(2 . l-2b) 
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✓ 



The resulting weighting or transfer function of the filter is given 
by : 

cp (s) 

T zx 



H(s) = 



f zz (s) 



+ f V 

v zx< s > 



(2.1-3) 



The use of this equation may be demonstrated by the example 
shown in Figure 1. Consider that x(t) has an autocorrelation 
function : 

1 - I T I 

cp xx (T) = - e 1 1 (2.1-4) 



The Laplace transform of the autocorrelation function is then 



^ (s) ■ 




(2.1-5) 



which yields the power density spectrum, 



W (0)) = — 
x 2 tt 



1 + (JO 



( 2 . 1 - 6 ) 



This signal can be shown to correspond to the random 
walk signal shown in Figure 2, where x(t) is a rectangular wave with 
values + ^ or - “ and with zero crossings located at event points 
which are Poisson distributed with an average frequency of crossings, 
per second. This time domain signal is not unique since many signals 

may have the same power density spectrum. 

Now consider a white noise source where, 

9 w (t) = 6(t) (2.1-7) 

where 6 (t) is the Dirac delta function, which has the following Laplace 
transform: 



<p (s) = 1 

vv 



( 2 . 1 - 8 ) 
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TIME 






FIGURE 2 

RANDOM WALK SIGNAL 
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The autocorrelation function of the signal plus noise is 



cp (s) = 9 (s) + 9 (s) + 9 (s) + 9 (s) 

zz T xx T w T vx T xv 



(2.1-9) 



Since the signal and noise are assumed to be uncorrelated the last 

two terms are zero and 9 (s) reduces to: 

z z 



C P Z2 (S) = V (s) +c f ) vv (s) 



1 - s 



+ 1 = 



2 - s' 
1 - s' 



( 2 . 1 - 10 ) 



The crosscorrelation function of the signal plus noise, and the signal 



is : 



9 (T) = 

zx 



(v(t) + x(t)) x(t + T)dt 



OC 

•i 



v(t) x(t +T)dt + / x(t) x(t +T)dt 



I 



= j x(t) x(t + T)dt = 9 xx ( t ) 



( 2 . 1 - 11 ) 



Since the noise and the signal are uncorrelated, 

1 



cp (s) = 
T zx 



1 - s 



( 2 . 1 - 12 ) 



Factoring Eq. 2.1-10 results in: 



9> (s) = 

T ZZ 



2 - S 2 = (V7 + s) ( V 2 - s) 

1 - S 2 (1 + s) (1 - s) 



(2.1-13) 



Therefore, from Eq. 2.1-1, 
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and, 



cp + (s) 
T zz 



V7 



+ s 



1 + s 



* 2Z (S) 



vr- 



1 - S 



Now, 



1 - s 



< Pz» (s) 1_ __ 

<P' 00 l - 6 2 (/T- S) (1 + s) ( V 2 - s ) 

6 ^ 



2.414 

Therefore, from Eq. 2.1-2a, 



^K <S) 






1 + s 



Cp (s) 
T zz 



2.414 1 + s 



(2.1-14) 



(2.1-15) 



Substituting into Eq. 2.1-3, 

1 



H(s) = 



2.414 



\1 + s. 



1 + s 



2.414 I \JT + s 



(2.1-16) 



The resulting filter weighting function, or impulse response, is: 



1 e- '/ Ft , t ^ 0 



H(t) =< 



2.414 

0 , t < 0 

This filter is physically realizable, since it is causal. 



(2.1-17) 
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2. THE BODE- SHANNON SOLUTION 



The frequency domain technique for the Wiener-Kolmogorov filter 
as developed by Bode and Shannon,^ is to convert the signal power 
spectral density to that of white noise, and then have a filter 
operate on this white noise . The development of this technique 
appears in several references 2,14,17 anc j discussed in Appendix 
B. The following results are obtained. 

Given W^((J0), the power spectral density of the signal x(t) ; 
and W^((Jo), the power spectral density of the noise v(t) , consider 
the transfer function of the following minimum phase (zeros and 
poles in the left half of the s-plane) filter: 



B(o)) = ^W x (cc) + W (tu) 



( 2 . 2 - 1 ) 



From a white noise input this filter will produce a power spectral 
density equal to that of the signal plus noise. As shown in Appendix 
B, the following transfer function can then be obtained for the non- 
causal optimal filter. 



h ((JO) = 

uo 



W (<u) B(u>) 

[W ((JO) +W (ao)] 

X V 



( 2 . 2 - 2 ) 



The preceding equations are represented in Figure 3. Since h^(u)) I s 
non-causal, let 



Ralph J. Schwarz and Bernard Friedland, Linear Systems 
(New York: McGraw-Hill Book Company, 1965), p. 334-45. 
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FIGURE 3 



BODE-SHANNON THEORY 



, t < 0 



8 cu (t) 



0 



V fc > 



, t * 0 



(2.2-3) 



where g^(t) represents the portion of h^(t) which exists for t ^ 0. 
This is equivalent to the spectral factorization of Eq. 2.1-1 and 
Eq. 2.1-2a. The desired optimum causal filter has the transfer 



function 

G (cd) 

H( "> ■ rw (2 - 2 - 4 > 



Now consider the example problem again. The signal, x(t), has 
a power spectral density 



W (cd) = 
x 



2tt 



1 



1 + CD 



2 



(2.2-5) 



and the uncorrelated measurement noise is white, with a power spectral 
density 

W v (w) (2.2-6) 



The power spectral density of the signal plus noise is therefore 



W (to) = W (tu) + W (uj) 





1 1 1 


= i_ 


r q 

2 + uj 2 


2TT 


2 + 1 

1 + CD 


2tt 


1 + <u 2 



(2.2-7) 



This may be factored 



W z (u)) = ^ 



( Vi~+ mo (/I- jou) 

(1 + p) (1 - jcu) 



( 2 . 2 - 8 ) 



Thus, the minimum phase transfer function required to produce this 
spectrum from white noise, using Eq. 2.2-1, is: 
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(2.2-9) 



B< ” ) - m 



V2~+ joj 
1 + jin 



The inverse of this transfer function will produce white noise from 
the signal plus noise* 

Now, substituting in Eq. 2.2-2, 



h (w) 

UK 7 



1_ I l ]( V2~ + joj j 1 
2n \ 1 + oj l\ 1 -f jin / \/2 tt 



2tt 



2 + in' 
1 + 0)' 



n/T + 



JSL 



(2 + in) (l + jin) 



'/tt 



2TT 



(VT- jin) (1 + jin) V2^ 



L— , ( _i + L. 

2.414 V2rr I \jz - jin 1 + j 






(2.2-10) 



By taking the inverse Fourier transform of Eq. 2.2-10, the impulse 
response, or weighting function, is then 




h (t) => 
in 



J_ 1 



I 

\/2n 



2.414 



STT 



2.414 



, t < 0 



-t 



, t S 0 



( 2 . 2 - 11 ) 



which is non-causal. The best causal impulse response, from Eq . 2.2-3, 
is 
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0 



, t < 0 



gw (t) = 



1 1 -t 

e 



>/2rr 2.414 



, t ^ 0 






The desired impulse response is therefore, 



and 



H(uu) = 



G (uj) 

UP 

B(cu) 



\ 



^2tt \ 2.414 I \ 1 + jw 



1 + jm 

JT+ P 



V 2tt 



2.414 JT+ p 



H(t) = , 



1 -V 2 t 

e 



2.414 



, t < 0 



, t ^ 0 



(2.2-12) 



(2.2-13) 



This filter impulse response is identical to Eq. 2.1-17, as expected. 



3. DISCRETE TIME SOLUTION OF THE 
WIENER- KOLMOGOROV FILTER 



The initial problem to be solved for digital computer simulation 
of the Wiener -Kolmogorov filter is the realization of the discrete 
signal, x(kT), based on the required power spectral density 



W (uu) = — 

X 2tt 



(2.3-1) 



H +® 

Since the power density spectrum of white noise is a constant, this 
power spectral density may be obtained by driving a linear system, 



G^(jw), by white noise, w(t) , so that 
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W (CD) = G.(ju>) G.(-juj) W (u>) = 
X 1 1 w 



1 



(2.3-2) 



2tt 

where 

G i<» - iVis 

It is extremely difficult, if not impossible, to obtain a 
discrete white noise source which is flat over all frequency. In 
practice, what needs to be done is to make the spectrum of the noise 
source flat over the bandwidth of The 3 decibel cutoff 

frequency for G^(jco) is one radian per second. The effective noise 
bandpass is given by 



1 + u> 



B = “ (1) = ”• radian/sec. = 1/4 Hz. (2.3-3) 

TT * 

where noise bandwidth is — times the 3 decibel cutoff frequency. 

The sampling rate for the discrete noise should be much greater than 
1/4 Hz., and 100 Hz. was chosen as convenient. The sampling period, 

T, therefore is 0.01 seconds. The spectrum of the noise should then 
be flat to 50 Hz. as shown in Figure 4. 

The autocorrelation function of this band limited white noise 
is then given by 



9 

1 T. 



. . 100 sin (100 ttt) 

, V / “ -I rid 'TTrr* 



(2.3-4) 



1 WW ' ' 100 ITT 

as shown in Figure 5. The variance of w(t), the noise sequence, is 
given by 



cp (0) = ct 2 = 100 = ^ 
T ww w T 



(2.3-5) 



Mischa Schwartz, Information Transmission, Modulation, and 
Noise (New York: McGraw-Hill Book Company, Inc., 1959), pp. 206-207. 
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FIGURE 4 

WHITE NOISE SPECTRUM 




FIGURE 5 

AUTOCORRELATION OF WHITE NOISE SOURCE 
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The white noise source is then simulated by a series of pulses of 
repetition period .01 seconds, whose amplitudes have a Gaussian 
distribution with zero mean and variance 100, followed by a zero 
order hold as shown in Figure 6. The recursive equation necessary 
for the computer simulation of the desired signal source is derived 
by taking the Z transform of the input-output transfer function of 
Figure 6. 



x (z) 7 


1 I 


1 - sT 

1 - e 


w(z) 


8+1 


s 




_\ 1 


1 J 




-T 




1 - e 








-T 





z - e 



Hence 

x(z) (z - e x ) = w(z) (1 - e T ) 



(2.3-6) 



and 

% 

zx(z) = x(z)e T + w(z) (1 - e T ) 



Therefore , 

x(k + 1) = e T x(k) + (1 - e T ) w(k); T = .01 sec. (2.3-7) 



To check that the signal does, in fact, have the desired power 

# 

spectral density, a subroutine entitled HARM was used to obtain the 
Fast Fourier Transform. The Fourier coefficients were then squared 
and plotted versus radian frequency, and the graph. Figure 7, obtained. 
Included on this graph is the desired power spectral density curve. 



International Business Machines Corporation, System /360 
Scientific Subroutine Package, (360A-CM-Q3X) Version II Programmer^ 
Manual, (White Plains, New York: International Business Machines 
Corporation, 1967 ) 
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FIGURE 6 
SIGNAL SOURCE 
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The program was run for 2048 samples and it can be seen that the 
power density spectrum approximates the desired spectrum. The 
computer program for generating the signal and checking its power 
density spectrum is included as Appendix C. 

Now that the proper signal has been obtained, a recursive 
equation for the filter is needed. The filtering system is drawn 
in Figure 8. Using the Z transform, the recursive filter equation 
is developed as follows: 



iizl 

z(z) 



1+/T 



s + -fZ I 



I 



1 


L -st| 


1 +/Z z 


r* ) 


1 / 


i - e- ^ T 


2 +VT" 1 z 


. e -^ T 



_L 

itfr 





(2.3-8) 



Hence 



x(z) 



z - e 



-VT 



T+\pr 



i-.-VT 



I 



z(z) (2.3-9) 



and 



z2(z) = &(z) e” ^ 2 T + - 1 



2 + 



V 2 



1 - e 



- \/2 T 



Z(k) 



(2.3-10) 



Therefore, the recursive equation is: 



!< 



k + 1) = e 



- \fl ' 



2(k) 



+ 2 +JT 



1 - e 



- V2” ■ 



Z(k) (2.3-11) 



The error is 

A 

E (k) = X(k) - X(k) 



(2.3-12) 
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DISCRETE WIENER -KOLMOGOROV FILTER 



The assumption is made that X(0) is zero. The sampling interval 
is 0,01 seconds. The measurement noise, v(t), is white, with a power 
spectral density, W^(w) = and is developed in an identical fashion 

to w(t), the driving noise for the signal generator. The two noise 
sources are uncorrelated. The mean square error and the variance 
of the error are then calculated, on a continuing basis, using the 
equations : 

1 n 

E (n) = ^ S E(k) (2.3-13) 

k=l 

E (n) 2 = — S (E (k) ) 2 (2.3-14) 

n k=l 

2 

ct = E 2 (k) - E (k) (2.3-15) 

e 

The computer run was for 30 seconds, i.e., 3000 samples, so that 
steady state is reached. The results are discussed in Chapter V, 

The computer program for the Wiener-Kolmogorov filter appears as 
Appendix D. 
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CHAPTER III 



BODE- SHANNON DISCRETE TIME FILTERING 



The Bode and Shannon frequency domain approach may be applied 
in the time domain directly, as developed by Schwarz and Friedland.^ 
The formulation is now in terms of summations, discrete signals, and 
discrete power spectral densities. This technique results in the 
following equation which is similar to the Wiener-Hopf integral 
equation, 

n 

^ v ( k > n ) = S h (n , j ) cp ( j ,k) (3.1a) 

zx , A zz 

J=0 

k - 0, lj 2 j .... n 



where 



c P^ v (k,n) = E z(k) x(n) 

ZX I 



(3.1b) 



is the discrete crosscorrelation function between z(k) and x(n), 
and 



<P 22 (J>k) - B 



z(j) 2 (k) 



(3.1c) 



is the discrete autocorrelation function of z(k). The unit impulse 
response of the filter is h(n,j). The notation h(n,j) denotes the 
response at time n due to an impulse applied at time j where j ^ n. 
Since the impulse response of the filter must be causal, this requires 
that h(n,j) equal zero for n less than j. Equation 3.1 holds for both 
stationary and non- stationary signals. When the signals are stationary 
then 

cp (k,n) = 9„ v ( n - k ) (3.2a) 

ZX Z X 

<P 2z (j> k ) = <P zZ < k - j) (3.2b) 
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h(n, j) = h(n - j) 



(3.2c) 



and Eq. 3.1 reduces to 



00 



9 (i) = E h(m) 9 (m - i) 

ZX - zz 

m=0 



(3. 2d) 



where 



i = n - k and m = n - j 



(3.2e) 



The power spectral density of a discrete signal, z(kT), is defined 



This is the discrete form of one of the Wiener-Khintchine equations. 
The prime denotes the power spectral density of a discrete signal. 

It should be noted that the z variable within the parenthesis denotes 
the Z-transform shifting variable. Proceeding as in Appendix B, a 
transfer function, B(z), is developed which will convert a signal with 



as : 



W z (z) = T Z cp^UT) 



00 




(3.3) 



a power spectral density of W^z) to a white noise process, u(kT) . 
This requires that 



B (z) B(z -1 ) w'(z) = 1 
z 



(3.4) 



Since B(z) is causal, its Z transform must be given by 



os 



B(z) = £ b(n) z 

n=0 



(3.5) 



Schwarz and Fried land 
of 1/T instead of T. 



12 

use this definition. Kuo uses a factor 
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If W^(z) can be factored into a product, 



w' (z) = W (z) W (z -1 ) 
z z z 



then 



(3.6) 



B(z) = (3.7) 

W z (z) 

For a white noise input 9^0 equals 6(j,k), and Eq. 3.1 reduces to 
/ 

cp (k,n) = h(n,j) for j ^ n 

ux 



cp^Ckjn) = 0 for j > n (3.8) 

The impulse response of an optimal causal filter for a white noise 
input, u(kT), and a desired output, X(kT), is then given by 



g(n, j) = cp (k,n) (3.9) 

U-X. 

For stationary processes Eq. 3.9 reduces to 

g(m) = cp (m) for m = n - k > 0 (3.10) 

ux 



Taking the Z transform of Eq. 3.10, with the aid of Eq. 3.3, yields 



g(z) 



W (z) 
ux 



(3.11) 



where W is the cross power spectral density of the white noise, u(kT), 
ux 

and the signal. The subscript _c denotes the causal part of the transform 
(terms involving z ^) . However 



i 



W 

ux 



(z) 



= W zx (z) B(z' 1 ) 



(3.12) 



so that 
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g(z) = 



(3.13) 



W^ x (z) BCz" 1 ) 



The optimum filter is thus 



H(z) = G(z) B (z) 



w' (z) B(z _1 ) 
zx 



B(z) 



(3.14) 



The procedure of taking "the causal part of" is analogous to the 
spectrum factorization done in Chapter II. 

Using the above technique, the previous example is now solved 
for discrete time signals. The signal and noise correlation functions 



are 



=K |k|T 



' XX 



cp (kT) = 6 (k) 
w 



(3.15) 



(3.16) 



The sampled power spectral density of the signal may be obtained by 
means of Z transforms and is as follows: 



W (x) = T Z 
z 



= T 



1 - kT 

2 6 



1 - e' 



(1 - ez b (1 - sz) 



(3.17) 



where e equals e . For the noise, the sampled power spectral density is 

w' (z) = T (3.18) 

v 

Therefore, the power spectral density of the signal plus noise input 



to the filter is 



w' (z) = w' (z) + w' (z) 
z x v 



= T 



2 - ez - ez 

(1 - Sz b (1 - ez) 



(3.19) 
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This may be written in the following manner for ease in factoring: 



w '( 2)= I£ (1 - bz'h (i - bz) 

z ' b (1 - ez 1) (1 - ez) 



(3.20) 



where 



b = 



i . / 



1 - e‘ 



(3.21) 



The transfer function needed to produce white noise from the signal 
and noise input to the filter is then, 

■1 



B (z) = V b/eT 



1 - ez 



v*> 



1 - bz 

Since the signal and the noise are uncorrelated. 



(3.22) 



and 



< f 2X (k) * < PK X <k) 



w' (z) - w' (z) 



ZX X 

Thus, from Eq. 3.13,L 



g(z) = 



w'(z) B(z _1 ) 

X 



1 - e' 



(1 - ez' ) (1 - ez) 



'/b/eT 



Simplifying and taking the causal part yields 

2 



g(z) = Vb/eT - 



1 - e 



be 1 - ez 



-1 



(1 - ez ) 

(1 - bz) 



J c 



(3.23) 



(3.24) 



(3.25) 



(3.26) 



The discrete time transfer function for the optimum filter is then 



H(z) = g(z) B (z) 



b_ 

Te 




(3.27) 
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The optimum filter diagram is shown in Figure 9. Since b, e, 
and T are constants, H(z) may be rewritten as: 



H(z) = C 



1 - bz 



-1 



(3.28) 



where 



C = 



b (1 - e ) 
Te (1 - be) 



(3.29) 



The recursive equation for the filter is developed as follows: 

A . 



H(z) = Sisl = 

n ^ z; z(z) 



1 - bz 



-1 



x(z) (1 - bz b = Cz(z) 



$(z) = bz ^ x(z) + Cz(z) 



x(k + 1) = b4(k) + CZ(k + 1) 



(3.30) 



The error is 



E (k) = X(k) - X(k) 



(3.31) 



For the specific example the coefficients for the filter, Eq. 3.30, 
become b = .869 and c = 12.45. These differ from the coefficients for 
the discrete Wiener-Kolmogorov filter, Eq. 2.3-11. One might expect 
that these results would be identical, but although Schwarz and Friedland 
consider both developments the discrepancy is not explained and remains 
a subject for further investigation. 
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CHAPTER IV 



THE KALMAN FILTER 
1. DISCRETE TIME SOLUTION 



The objective in this chapter is to examine the filtering 
problem from the state point of view as developed by Kalman . ^ ^ 5 ^ 
In the previous cases, the signal was characterized by either the 
autocorrelation function or the power spectral density. Here the 
signal is characterized by the output of a linear dynamic system 
driven by white noise similar to Bode- Shannon . This presents no 
problem, for if the power spectral density is rational, it can 
always be represented as the output of a linear system driven by 
white noise. As in Chapter II, the linear system has the transfer 
function, 



G 1 (s) = 



1 

s + 1 



(4.1-1) 



Since this is a first order system, the observation matrix is unity 
and scalar notation can be used. The corresponding differential 
equation is 

= -x(t) +w(t) (4.1-2) 



Then, from state variable theory, 

^ = e ^ and F = 1 - e 



and the differential equation may be expressed in equivalent recursive 
form: 

X(k + 1) = e' T X(k) + (1 - e“ T )W(k) (4.1-3) 
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The signal model is outlined in diagram form in Figure LO; w(t) 
and v(t) are independent gaussian random processes with zero means 
and covariances given by: 

cov [w(t), w(t)] = Q6 (t - t) 
cov [v(t), v(t )] = R6 (t - t) 
cov [w(t) , v(t )] = 0 

Q and R are the variances of the signal white noise source and the 
measurement noise, which are set equal to unity in order to duplicate 
the model for the Wiener filter. The noise sources are simulated in 
the exact manner as described in Chapter II, Section 4. 

The Kalman filter^ is shown in Figure 11. G(k) is an adjustable 
gain which minimizes the mean square error, 



2 1 n 2 
E(n)=~ S (E(k)r 

n k=l 


(4.1-4) 


where 




E (k) = X(k) - X(k) 




G(k) is determined by the following equations: 




G(k) _ mPLZ-11 

P(k/k - 1) + R 


(4.1-5) 


P(k/k) = (1 - G(k) ) P(k/k - 1) 


(4.1-6) 


P(k + 1/k) = P(k/k) + QD 


(4.1-7) 


where P(k/k - 1) denotes P at time k given the value at 


t ime k - 1 . 
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KALMAN SIGNAL MODEL 
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P (k/k- 1) = E 




l 

k 



k-1 

E 

n=0 



Z(n+1) - Z(n+l/n)J 



P(k/k) = E 



Z(k) - Z(k/k) 



2 




Z(n) 



Z (n/n) 



2 



also, 

QD = r 2 Q 



It is necessary to define P(l/0) in order to start the recursive 
process , 

It is important to note that once P(l/0) is specified the gain 
equations, Eqs . 4.1-5, 4.1-6, and 4.1-7, are recursive and do not 
depend on the observations Z(kT). They depend upon the variance 
of the signal driving noise, the variance of the measurement noise, 
and the characteristics of the linear system used for signal generation. 

From Figure lj, the recursive equation for the Kalman filter may 
be expressed as 



X(k) = §X(k - 1) + G(k) 



Z (k) - §X(k - 1) 



(4.1-8) 



The error is 

E (k) = X(k) - X(k) 



(4.1-9) 



As in the Wiener simulation, the sampling interval is 0.01 
seconds, $(0) is the initial estimate at time t = 0 and the length 
of the computer run is 3,000 samples. The computer program for 
the Kalman filter simulation is included as Appendix E. 



44 



2. EQUIVALENCE OF KALMAN AND WIENER- KOLMOGOROV FILTERS 

IN THE STEADY STATE 

It has been shown by C, E. Hutchinson^ that the continuous 

Kalman and Wiener filters are equivalent in the scalar case after 

the steady state condition is reached, i.e., G(k) is a constant. 

This equivalence is developed as follows. 

a 2 

Given a message spectrum, W (U)) = ■ ■ 

x bV + 1 

2 

a noise spectrum, W^(tt)) = c 

and a crosscorrelation function, 9 (t ) = 0 

T xv 



The optimal causal filter, in the Wiener sense, is then 



where 



H(s) = (k 2 - k L ) 



1 

s + k 2 



\ = 1/b 

k 2 = l/bc ^a 2 + c 2 



(4.2-1) 



Now consider the Kalman filter. The signal model is, in general, 
defined by: 

x(t) = Fx(t) + w(t) (4.2-2) 

z(t) = Mx(t) + v(t) (4.2-3) 

The best estimate of x, designated :x, is obtained from the continuous 
filter equations which correspond to the discrete equations 4.1-5, 
4.1-6, 4.1-7, and 4.1-8. 
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x(t) = Fx(t) + g(t) [z (t ) - Mx(t)] 



(4.2-4) 



where 

g (t) = p(t) m t r _1 

P(t) = FP(t) + Pf T (t) - P(t)M T R _1 MP(t) + Q 

and 

cov [w(t), w (T ) ] = Q6(t - t) 
cov [v(t), v(t)] = R6(t - t) 
cov [w(t) , v(t)] = 0 



If a one dimensional problem is considered with 



F = -1/b 




2 



the Kalman solution becomes 



x(t) = -1/b x ( t ) + g(t) [ z (t ) - x ( t ) ] (4.2-5) 

where 

g(t) = P(t)/c 2 

P(t) = - (2/b)P(t) - (l/c 2 )P 2 (t) + a 2 b 2 



In the steady state, when P(t) is identically zero, the equation for the 
Kalman filter reduces to 

x(t) = -1/b x(t) + 



-1/b + 1/bc 




+ c 



z(t) 



X(t) 



(4.2-6) 
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The transfer function relating the estimated signal, X(s), to 
the signal plus the noise, Z(s), can be expressed as 



ksi . 

Z(s) 


■ (k 2 ' V s A 2 < 4 - 2 - 7 > 



which is identical to Eq. 4.2-1. This result can be used as a check 
on computer programs for the Kalman and Wiener -Kolmogorov filters. 
When G(k) has reached a constant value, the mean square errors should 
be identical. 
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CHAPTER V 



RESULTS AND COMPARISONS 

The continuous Kalman and Wiener-Kolmogorov filters have been 
shown, in Chapter IV, to be equivalent in the steady state condition, 
that is, when the Kalman gain is constant. This result is verified 
experimentally, as shown by the computer curves of Figures 16, 20, 24, 
and 28. The steady state gain of the Kalman filter reaches a value 
of .0041335 as expected. As soon as the gain reaches its steady 
value, the curves of mean square error versus time (Figures 18, 19, 

22, 23, 26, 27, 30, and 32) are identical to each other and to the 
Wiener-Kolmogorov filter curves (Figures 14 and 15). 

In order to start the recursive process, the Kalman filter 
requires a priori knowledge of the initial error covariance, P-^/q* 

Four different values of P^q (0*0, .414065, 1.0, and 10.0) were 
chosen to compare the effectiveness of the Kalman filter for different 
values of this term. The value of P^q = .415065 was chosen because 
it gives the steady state gain (.0041335) at the start of the computer 
run. With this P^q t ^ e performance of the Kalman and Wiener-Kolmogorov 
filters are identical (Figures 13-15 and 21-23). 

In simulating the non-s tationary signal, two different initial 
values were chosen; one to correspond to a small initial error in 
estimation and the other to correspond to a large error in the initial 
estimation. In a tracking problem, the first case would represent 
the situation where a target is picked up at long range. The second 
case represents the situation where a target suddenly appears at 
relatively short range. For small initial errors, the Wiener-Kolmogorov 



48 



and the Kalman filters performed essentially the same (Figures 13, 
14, 17, 18, 21, 22, 25, 26, 29, and 30). For large initial errors 
the Kalman filter, with large initial covariance of error, P^q, 
gives the best transient response (Fig. 31). For large initial 
error, if ec L ua l to zero, the transient response of the 

Kalman filter takes somewhat longer time to settle than the Wiener- 
Kolmogorov filter as shown in Fig. 19. 
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CHAPTER VI 



CONCLUSIONS 

The responses of the discrete Kalman filter for various initial 
error covariance and the discrete Wiener -Kolmogorov filter have been 
compared for a non-stationary random input with defined statistical 
properties. It was necessary to develop a computer program for the 
simulation of the desired spectrum from a white noise source. The 
discrete time Wiener-Kolmogorov filter was also programmed for the 
digital computer. 

It has been shown that when the gain of the Kalman filter is 
equal to its steady state value the performance of the Wiener- 
Kolmogorov and Kalman filters are identical. For large initial error 
the Kalman filter (with a large initial gain) shows the best transient 
response. When the initial error is small the Wiener-Kolmogorov and 
the Kalman filters respond in as essentially equivalent fashion. 
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APPENDIX A 



DEVELOPMENT OF THE WIENER-HOPF INTEGRAL EQUATION 



The Wiener-Hopf equation gives the relationship between the 
filter weighting function, the autocorrelation function of the 
input to the filter, and the crosscorelation function of the filter 
input with the desired filter output. The problem is described in 
Figure 1; x(t) is some signal, v(t) is white noise, x(t) is the 
estimate of the signal, and e(t) is the error in estimation; reduce 
e(t) to a minimum in the mean squared sense by a correct choice of 
h(t), the impulse response of the filter. The error by definition is 



e(t) = x ( t ) - x(t) 



Squaring both sides yields, 



e 2 (t) = x 2 (t) - 2x(t) x(t) + x 2 (t) 



(A. 1) 



(A. 2) 



The output, x(t), and the impulse response of the filter, H(t), may 
be related by means of the convolution integral. 



x(t) = 








00 




00 


A 2 . v 

X (t) = 


1 H(t) z(t - T)dr 




J H(T 1 ) Z(t - T 1 )dT 1 




J 

L. -00 J 




✓ 

- -03 



1(t) x(t - t) dT (A. 3) 

where T is the variable of integration. Squaring Eq . A. 3 yields 

(A. 4) 

L L ! 

Substituting Eq. A. 4 into Eq . A. 2 gives, 

— 00 -» 

e 2 (t) = x 2 (t) - 2 J H(t) z ( t - T)dT x ( t ) + 

_ -00 

(A. 5) 



p 00 






f H(T) z(t - T)dT 
^ -00 




HCt^) z(t - T 1 )dT 1 
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By definition, the mean square error is, 

T 

e 2 (t) = lim 



1 _ f 

1 im 2 'j] I 

J T 



e (t )dt 



(A, 6) 



Substituting Eq. A. 5 into Eq . A. 6 yields: 



e 2 (t) 



lim ^ x 2 (t)dt 
T^co 2 T ' 



2 lim ^ 
T->oo z 1 



1 o 

j ic j 

-T •'-e 



H(T)z(t - T )x ( t )dT 



+ ^ 2 T 1 dt 

-T 



T o 

i> I 



H(T)z(t - T )dt 



H(T 1 )z(t 



T l )dT l 



(A. 7) 



The definition of the input autocorrelation function is 

T 



cp (t) = lim 
zz T_»<n 



•i f 

2 T J 



z(t)z(t + T)dt 



(A. 8) 



Similarly, the crosscorrelation function between the input and the 
output is 



A - . I 

2 T 



9 Z ^(T) = m l.irn ^ / z(t)x(t + t) dt 



(A. 9) 



-T 



and the autocorrelation function of the output is. 



Cp xx (T) “ ii 4 m 2 T r + T > dt 



(A. 10) 



Now, making use of the correlation functions, and interchanging the 

order of integration in Eq. A. 7 yields, 

. 00 

e*(t) = cp (0) -2 f H(T) cp (T)dT + 

XX I zx 



co co 

j H(T )dT 
-00 



1 

J H(t 1 ) cp 



zz< T ' T l )dT l 



(A. 11) 
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To determine the system function that minimizes the mean square 
error, assume that a solution exists and denote this solution as 
H m (t). Now construct a weighting function: 



H(t) = H (t) + 6H, (t) (A. 12) 

m o 

where H^(t) is the assumed solution, Hg(t) is any arbitrary realizable 

weighting function, and 6 is a parameter which may be varied to test 

whether H (t) is the solution. The construction is such that at 6 
m 

equal to zero, the derivative of mean square error with respect to 6 

must be equal to zero. By setting this derivative equal to zero for 

6 equal to zero, a condition for H (t) which must be satisfied in 

m 

order for it to be a solution, is specified. 

Substituting H(t) as given by Eq . A. 12 into Eq. A. 11 and differenti- 
ating with respect to 6 yields: 







d 


? 

Le (t)J 


d6 



2 ^ H 6 (T)tp zx (T)dT + J V T)dT / WW 1 • V dT ! 



00 



00 



+ J H fi (T)dT 



I H (T ) CP (T - T )dT 
j m 1 zz 11 

-00 



+ 26 J H 6 (T)dT J H 6 (t) CP zz (T - T l )dT 1 



(A. 13) 



Because of the even property of autocorrelation functions of stationary 
signals , 

• T) ■ ^ 2 <t • V <A - 14) 

which leads to 
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/v T)dT / VV f S2 (T - T i )dT i 



00 00 



V T)dT H m (T l )c P ZZ (T - T l )dT l 



(A. 15) 



Substituting Eq . A. 15 in Eq, A. 13 and setting 



d _L _e lt)j _ 

d6 



= 0 at 6 = 0 



yields : 



2 H 6 (T)dT 



H (T )cp (r - T )dT - CP ( T ) 
m 1 zz 1 L zx 



= 0 (A. 16) 



Since Hg (T ) is a realizable weighting function, it must be zero for 

values of t less than zero. The only way that Eq . A. 16 can be 

satisfied for T ^ 0 is that the factor in the brackets be equal to 

zero. Thus, 

00 

( H (T ) cp (t - T )dr. = cp (T), for T = 0 (A. 17) 

J m l zz \ 1 T zx 

-CO 



Equation A. 17 is the famous Wiener-Hopf integral equation. 

Strictly speaking, the H^Ct) that satisfies Eq. A. 17 produces a 
stationary value of the mean square error and does not necessarily 
ensure a minimum. It can be shown that the solution of Eq. A. 17 does 
yield a minimum, by considering the second derivative of the mean 
square error with respect to 6. Differentiating both sides of Eq. A. 13 
yields 




00 



00 



2 J V T)dT / W V (T ‘ V dT i 
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This is the same as twice the mean square value of the noisy input 
signal filtered by a weighting function, H^(t). Since this mean 
square signal can never be negative, this shows that the second 
derivative of the mean square error is always positive, hence the 
solution of Eq. A . 17 corresponds to a minimum. 
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APPENDIX B 



DEVELOPMENT OF THE BODE- SHANNON SOLUTION 



Bode and Shannon used frequency domain techniques to solve the 
filter problem by converting the actual signal power spectral density 
to that of white noise, and then operating on this white noise signal. 
The solution may be described by considering Figure 1; x(t) is the 
signal, v(t) is white noise, x(t) is the estimate of the signal, and 
e(t) is the error. The assumptions are that, (1) input signal and 
noise are uncorrelated and stationary, and (2) all time functions are 
Fourier transformable. 

The impulse response of the filter, h(t), is chosen on the basis 
of minimum mean square error. The error signal is defined as 



e (t) = x(t) - x(t) 
and its Fourier transform is, 

E (u>) = X(UQ) - X(to) 

= ^X(ou) + V(u))j h(u)> - X(cjo) 

= [h(w) - lj X(w) + h(uu)V(aa) 



(B.l) 



(B.2) 



Since the signal and noise are uncorrelated, the power spectral density 
of the error is 



W (uo) 
e 



1 

2tt 



|h(cu) - 1 | 2 W x (u>) + |h(eu) | 2 W v (w) 



(B . 3) 



Hence the mean square error is 



e 2 (t) - ± 



|h(u>) - l| 2 W x 0») + | h (a)) | 2 W v (w) 



dtu (B .4) 



This is to be minimized by the proper choice of h(iu). 
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Note that to obtain a solution, the power spectral densities 



W (u>) and W (w) must be specified. Let 
x v 



h(to) = 

= C(to) (CosQ(to) + j sin0(to)) 
Substituting into Eq. B.l and simplifying yields, 



(B.5) 



e 2 (t) 



-i- f 

2tt J 



00 

/ <|c 2 (to)W (to) + 


2 ll 

C (to) + 1 - 2C(to)cos (9(to) ) 


l ( V 


J J 



' W x (cu)du} 



(B.6) 



Since C(to), W (to), and W (to) are nonnegative this expression is 

X V 

minimized when cos 0(u)) has its largest value; that is, when 0(u)) 
equals zero. Thus, 



2 



e 



opt 



(t) 



1_ 

2tt 



00 



/ 


C 2 (t0) 


W (to) + w (to) 1 

V x J 


- 2C (to )W (to) + W (to) 


-*oo 






_ 



du) 



(B . 7) 



In order to find the C(to) that minimizes this expression, complete 
the square of the terms in the bracket. Thus, 



2 

e opt 



(t) 



CO 



1_ 

2tt 




:(to) \/w~ 



(to) + W x (to) 



v*> \ 2 

Vw v (t0) +W x (tO) J 



W (to)w (to) 

X V 

W (o)) + W (0)) 

X v 



Eq. B.8 is minimized when 

W (u)) 

C(to) = 2 

W (to) + W (to) 

X V 



(B.8) 



(B.9) 
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with minimum mean square error given by 

CO 



e . (t) 

min 




W (uu)W (uu) 
x y 

W (uu) + W (uu) 

X V 



du) 



(B.10) 



Since 9(uu) equals zero, 
h(uu) = C(uu) * — 



W (uu) 
x 



W (uu) + W (uu) 
x v 



(B . 11) 



The impulse response is therefore given by 

CO 



h(t) = — 



W (w>) 
x 



2tt I W k ( “> + V“> 



juut , 
e J duu 



(B.12) 



which, in general, does not vanish for t < 0. Thus Eq „ B.ll is 
generally not physically realizable. 

The minimum phase filter which converts a white noise input to a 
spectral density equal to that of the signal plus noise, must have 
a transfer function of magnitude 

B(u )) = V W (u>) + W (oj) 

X v 



Since Eq . B.ll gives the best noncausal operation on the signal plus 
noise, the best noncausal operation on white noise inputs requires that, 






W (CJU) y ; W 0 ), ) + W (0)) 
X v X v 

W (ao) + W (ao) 

X V 



(B . 13) 



This however, is still noncausal. The causal filter is then constructed 
as follows: 



, t < 0 



8«. (t) = < 



^^(t) , t * 0 



(B.14) 
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The desired causal filter for operation on signal plus noise 
inputs then has a transfer function given by: 



H «"> = g u> ( uj) 



1 

V w (cu) + W («)) 

X V 



(B.15) 
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APPENDIX C 



POWER SPECTRAL DENSITY PROGRAM 
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FORTRAN IV G LEVEL 0, MOO 0 



POWER SPECTRAL DENSITY PROGRAM 
MAIN OATE 



67328 



18 / 15/48 



0001 



OC C2 
o 0 0 3 
<UC4 
000 5 
CCC6 
OCC 7 

ooc e 

0Co9 
CO 10 
00 11 
oO 1 2 

0013 

0014 

0015 

0016 
0017 
GO 18 
Ct 1 9 
)C2C 
0021 
DC 2 2 
UC2 3 
0( 24 
CC 2 5 
OC 26 

0027 

0028 
0029 
003C 
00 3 I 
LL 3 2 

0033 

0034 
OC 3 5 
OC 36 
0 0 37 

0038 

0039 

0040 
0C41 
0042 
Co 4 3 
0^4 4 

0045 
C C 4 6 

0047 

0048 
OC 4 9 

0050 

0051 

0052 

0053 

0054 

0055 

0056 
CC 57 



c 7=iRlf A biflE^l<i}k E A£ 2 Input ar r a y - to - har m-m u s t 

C BE A POWER OF TWO 

CCMPLEX*0 A<2048,1,1> 

DIMENSION SI 10241 , INVU024I 

DIMENSION M< 31 ,X< 41001 tB<20aWBX< 200) ,TIMEX! 2001 
K K= 1 1 
T=C • 0 1 

PIT=6. 2831053 

00=1. 0/T * 

D E V= SQR T I QC I 

R=1.0/T 

DEX=SCRT!R) 

PH I = E XP ( -T ) 

GA MMA= 1.0-PHI 
I X=4C97 3 4 7 51 

CALL GAUSS! IX,DEV,0.0,W) 

XIl)=GAMMA*k — - — 

MI 1 ) =KK 
M I 2 ) =0 
MI 3)=0 
N 1=2**KK 
N 2= 1 



N 3= 1 
CC 2 
DC 2 
OC 2 



I 1 = 1 ,N1 
12=1, N? 

13=1, N3 

2 A( I 1, I2,I3) = I0.C,0.0) 

SIGNAL IS NOW GENERATED 

DC 3 I 1=1 ,N1 
DO 3 12=1, N? 

Of 3 13=1, N3 

CALL GAUSS! IX, DEV, 0.0, W) 

XII 1+1)=PH I*X< I 1)+W*GAMMA - — 

3 A ( Il,I2,m = X< ID 

CALCULATE Mp AN AND MEAN SQUARE OF SIGNAL 
AX=0 . 0 
B Z = C . 0 

DC 50 1=1, M 
5C AX=AX+X( I ) 

XB AR = AX/M 

DO 51 1=1, N1 

5 1 BZ=BZ+< XI I )-XBAR)**2 
XMS=BZ/N1 
WPITF16,92)XEAR 

92 F0RMATIT6, f SIGNAL AVERAGE VALUE = «*T32,Ei3.6) 

WRITEI6,95)XMS 

95 FORMAT! T6, *MFAN SQUARE SIGNAL = »,T30,E13.6) 

OBTAIN THE FAST FOURIER TRANSFORM 

CALL HARM1 A , M , 4 N V * 

N X= 3 5 

NORMALIZE AND PREPARE INPUTS FOR PLOTTING 
DC 42 1=1, NX 
B I I ) =C A8S ( A I 1,1,1) >**2 
B X I I ) = ( 61 I )*2.G) /I N1*QD) 

TIMEX! I )=( I I-1)*PIT)/!N1*T) 

WRITE(6,63)TIMEX< I ),BX< I ) 

63 FORMAT! T6> r E 13^6, T2C , E 1 3, 61 - — 

42 CONTINUE 

READ <5, 105)1 I TITLED), I = 1,12) 

105 FORMAT { 6A0 ) 

CALL DRAW! NX, TIMEX, BX, 0,0, * • ,ITITLF, 0.0,0. 0,0,0,0?0, 8,9, I, LA) 

END 
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APPENDIX D 



WIENER-KOLMOGOROV PROGRAM 
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fCklRAN IV C LEVEL 0, MOD C 

m\ 

CCC 2 
COCA 
OCC 5 
C0C6 
0007 
COOP 
COLS 
001C 
CC 1 1 
00 12 
OC 13 
00 1 A 
CC 1 ^ 

OC u 
lO 1 7 
OGie 
C01S 
002C 
0021 
0022 
0022 

^C 2 A 
0025 
CC 26 
C 0 2 7 



WIENER-XOLMOGOROV FILTER PROGRAM 

MAIN DATE = 67330 



14/00/57 



CO 28 
CC2S 
CC3C 
GO 3 1 
CO 32 
CC 3 3 
CC3A 

CC 3 5 
CO 3 6 
CC ? 7 
00 3 8 
CC ?s 
C 0 AC 
C041 
CC A2 
A 3 
C C A A 
CCA 1 ' 
OCA 6 
C > 7 
C 3 A 6 
CCAS 

CC 50 
C 0 5 1 
CC 5 2 
0^5 3 
005A 
L- o 5 5 
CL 56 
0 0 5 7 
CC58 
C "5S 
r O 6C 
CO 6 1 

c 0 6 2 
GC 6 3 



BlffettiMWIASL > XMAT I 310 0 1 i6 < 3 1 QG TWO 3-1001 vT 1M€ I 3100I 

D I MENS 1 CN EMS( 31C0) tYXX( 310) • ES< 3100 I 

N = 3 C C C 
NA=N+5 

I X=A CS 7 3 A 7 5 1 — — 

T=0 • C 1 

QC= 1 • C/T 

CEV=SCRT!GC) 

R=1*0/T — - 
OE X= SCR T ( R > 

EMS ( l)=C.O 
PH I = E XP ( — T ) 

GAMMA=1.0-PH1 

RCT=SCRT!2.C) 

XPC=RCT *T 
PHK= EX P (- X PC ) 

GAM8=1*0-PRK — - - 

MR I TE ( 6 , 20 ) 

20 FORMAT (T6, *K*,T16,*XlK)» , T36 , • XH A T < K ) ' ,T58, *Z! K > * , T80, * EMS ! K ) * ) 
Xh A T ( 1>=0.C 

CALL GAUSS! IX, DEV, C.O,Wl 
X( 1 )=GAMMA*M 
C GENERATE SIGNAL 

DC 10 K= 1 t N A 

X! K + 1 ) =PH I *X ( K ) + M*G AMM A 

10 CONTINUE 

C ACD MEASUREMENT NOISE AND 

C COMPUTE ESTIMATE Of SIGNAL 

DC 11 K= 1 , N A 

CALL GAUSS! IX, CEX, 0.0, V) 

7(K)=X!K)+V 

XRAJ4KM ) = PKK*XH4Tm46AW8*444C|^n^<LU2.O*R0Tl ) 

E<K>=XHAT!K)-X!K> 

11 CONTINUE 
KT=2 

C COMPUTE MEAN AND MEAN SQUARE ERROR 

DO 85 K= 1 , N 

Y A = C • 0 

ve*c,o 

DC 25 1 = 1, KT 
YA = Y A + F ( I ) 

25 CONTINUE 
YPAR = YA/K T 
DC 26 1=1, KT 

Y 6= Y P+ ( E ( I )- YBAR )* + 2 

26 CONTINUE 
EMSCK+l ) = Y6/KT 

KT— K T+ 1 — 

ViPITE!6,30)K,X(K)*XHAT(K),ZCK) , EMS ( K ) 

30 E0RMAT(2X,IA,5X,Ei3.6,9X,E13.6,9X,E13.6,9X,El3.6) 

85 CCNTINOF 

C PREPARE ITEMS FOR PLOTTING 

CC 7 C 1=1, N 
TIME! I ) = l-l 
7C CONTINUE 

00- -1 10 -I=L r * — — 

110 E S ( I )= E ( I M*2 
K C= 1 
NPT= 3Q0 

DC 43 I=1,NA,10 
YXX( KQ) = EMS( I ) 

A3 KC = KQM 

REAC(5,1C1 ) ( ITITLE! I ) ,1 = 1,12) 

101 FCRVATI6A8) 



CAIL CRAM! NPT ,T1ME,YXX,0,C, * 

END 



* , ITITLE, 0*0, 0*0,0, 0,0, 0,8, 9, 1 , L A ) 



APPENDIX E 



KALMAN PROGRAM 
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KALMAN- PIJ^TBR PROGRAM 



FORTRAN IV-4 LiVH 0^ MOD 0 MAIM OATF * 6 733L 16/20/48 



m 

OOC 3 
OOC A 
OOC 5 
OCC 6 
QCC 7 . 






BHtSfl HlltilH) 


OI?ENSION‘illi6o) .XHAT(3100»,EOI00I,Z( 3100) ,T I ME < 3100 I , EMS ( 3 100 I 
DIMENSION YXX< 310 ) ,GX( 3101 »PK?3 1001 , PI 31*01 ,G( 31001 . ES ( 31 00 » 


T=C.01 ~ 

N= 3000 

NA = NU5 . _ _ ... _ . _ 


OCC 8 






Q0= 1 *0/T 


OCC 9 






DEV* SCR T ( Q C ) 


0010 






ft~V.O/T 


OCll 






DE X = SORT { R ) 


OC 1 2 






I X = AC973A7 51 


OC 1 3 






EMS( 11 = 0*0 


COM 






PHI=EXP(-T ) 


0015 






GAMMA=1. O-PHI 


C 0 1 O 






XHAT ( 11*0. C 


0 C 1 7 






q-gamma*gamma*<dev**2 ) 


0018 




— 


pi H = io*n * 


C C 1 9 






N P = N A + 2 




c 




generate GAIN secuence 


CC 20 






CO 57 K= 1 t NR 


0^21 






G(K)=P(K)/(P(K|4-R) 


G 022 






PK ( K 1 * ( l.C-GtK) ) * P ( K ) _ . 


C rt 2? 




57 


P(K + n*PHI*PhI*PK(KUQ 


•^?8 






UR IT E ( 6 , 2C I 


uG?5 




20 


FORMAT! T6, »K » ,T13, * XT K ) • , T 3?, • XHATL K ) ' ,T50, • 7< K ) * , T70 r • EMS ( K 1 * ,T90 








1 , f G ( K ) 1 ) 


0026 






CALL GAUSS I I X ,CEV , 0.0 , U ) 


OC 2 7 






X tTT= GAMMA * V “ ~ ~ 




c 




GENERATE SIGNAL 


0028 







£C 1C K=1,NA 


OC 29 






CALL GAUSS! 1 X,DEV,O.C,W) 


OC?C 






XIK+11 = PHI*X(K )+W*GAMMA 


oon 


- 


10 


CONTINUE 


0 32 






CALL GAUSS! 1 X .DE X f 0. 0 , V ) 


JO 3 3 






zm=x< n + v 




c 




AXD MEASUREMENT NOISE 




c 




CCMPUTF ESTIMATE OF SIGNAL 


GcB'i 






DC 11 K=1,NA _ . - 


(■ 3 9 






CALL GAUSS ( 1 X,CfX t C.O,V) 


0. 36 






Z(K+I )=X(K+l)+V 


0037 






XHAT!K4-H = PHI*XHAT(K)+G<K + 1>*I7IK4-1 )-PHl *XHAT( K ) I 


0*138 






E(K)=XHAT(K)-X{K) 


OC 39 




1 1 


CCNTINUE 


o:ao 






KT=2 




c 




CALCULATE ME AM AND MEAN SQUARE OF ERROR 


LC'il 






CC 85 l< = ) » N 


C : 8 2 






Y A = tj # 0 


o; a? 






YR= : .C 


OC A A 






CC 25 1=1. KT 


CO A 5 






Y A = Y A ♦ E { 1 ) 


OCA 6 




25 


CONTINUE 


0 C A 7 






YPAR=YA/KT 


OC A 6 






DO 26 1 = 1. KT 


00 A 9 






Y H = Y fi + ( E< I )-Y8AR)**2 


00 50 




26 


CONTINUE 


C"M 






FMSI Kf 1 I = Y P / K 1 


Oi. 52 






K T=K T* 1 


C.5 3 






WR ITE ( 6 , 2" )K ,X(K 1 , XHATIK ) , Z(K I ,EM$< K) ,G(K) 


0 '54 




30 


FORMAT ( I4»5X,E13.6,9X,E13*6,9X,E13»6,9X,E13 •6.9X.E13.6I 


c : ^ 6 




85 


CONTINUE “ ^ - — 




c 


PRFPARF ITEMS FCR PLOTTING 


C 5 






DO 1C 1=1. N 


^,57 






T I ME ( I 1 = 1-1 


^ . c p 




70 


CONTINUE 


Ot 59 






00 110 1=1, N 


l. 0 6 v. 




1 10 


ESm=EU 1**2 


006 I 






KQ-1 


00 6 2 






N* P T = 3 0 C 


0063 






DO 4 3 1 = 1, NA ,1C 


>„',6A 






YXX1.KQ ) =EM $ ( I) 


0^,6 5 




A3 


k‘Q=kQ+ 1 


>;V^ 






READ! 5, 1C 1 )( IT ITl E (I ) ,1 = 1 ,12) 
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*f 





KALMAN FILTER PROGRAM 


FOWTUAN 


IU C- L £ V 6 L— 0-»— MOO 0 - - -MALN- - DATE »67*M - M>/g0/48 


8:U 
0 0 6^ 
c J7c 
CO 71 
C )7? 
0072 
0)74 
CO 7 5 


101 ^[t^fAS6?dpT # TiM€ f VXXtOt©^r^ ! T ITL E , 0 . 0 ,0 .0,0 ,0,0,0, 8,9, 1, LA) 

c r’44 I = 1 , N A , 1 c 

GX ( KL ) -G ( I 1 

44 Kl =Kl* 1 

ORAwkp'f "TTKTITLE, 0.0,0. 0»0,0,0,0,8,9,i,LA) 

END 



88 



APPENDIX F 



NO FILTER PROGRAM 



89 



NO FILTER PROGRAM 



FCPIFAK IV G LEVEL C, *<00 0 

88U 

OOt 3 
OCC4 
OCC 5 
oorfc 
0017 

core 
GOES 
CGIC 
0011 
00 12 
CC 1 3 

0014 

0015 
OC 1 6 
0017 



MA IN 



DATE = 67330 



1 3/54/55 



1 xVliiSI, XH ATf3rl£G \ t ttH 3100 ) ,T IM€ ( 3iOO) 

DIMENSION EMS! 3 ICO) ,YXX! 310 1 , W! 3 100 1 , V ! 3 100 1 t E S( 3100 ) 



N=3000 
NA=N*5 
T = 0 • G 1 
QC = 1 • 0 /T 
DEV=SQRT!CC) 
R=1.G/T 
OEX^SS«T!R 1 
IX=4C9734751 

x< n=o.o 
>=c. 



00 IP 

001 <5 



00 2C 
0021 
00 2 2 
0023 
00 24 
00 2 5 

0026 
OC 2 7 
0C2P 
0C2S 
00 3 C 
OC 3 1 
0032 
00 3 3 
CO 34 
OC35 

0036 

0037 
00 3 P 
Cl?S 
00 4 1 
OC 4 1 

0042 
CO 4 ? 
00 44 
0 C 4 5 

0046 

0047 

0048 
OC 4 S 
C 0 5 C 
00 5 1 
00 5 2 
Ou53 
0054 
00*5 
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1C 



ems ( 1 i=c.o 

PH I=EXP ( -T 1 
GAMMA=1 .C-PH I 
WRITEI6.2C ) 

20 FORMAT! T6, • K • , T 1 4 , • X { K ) • , T34 , • Z ! K ) • , T 54 , • EMS (K 1 • f T 74 f • W ( K 1 • , T94, 
lr»VTK M \ 

CG 11 K= 1 , N A 

CALL GAUSS(IK.OEV,C.O,W(K) 1 
GENERATE SIGNAL 
ACD MEASUREMENT NOISE AND 
COMPUTE ESTIMATE OF SIGNAL 
Cf 10 K - 1 * N 

x < k+ n = PHi*xiK i + mk+d ♦gamma 

£ ALU G4LlSSUXr&€^rO^**LKn 
Z ( K 1 = V ( K 1 ♦ X ( K 1 
E(K1— Z(KI — X(KI 
CCNTINUF 

COMPUTE MEAN AND MEAN SQUARE ERROR 
KT=2 
CC 85 K= 1 t N 

Y A = 0 • C 



CC 25 1 = 1 » K T 

Y A= Y A *E ( I ) 

25 CCNTINUE 
YBA«=YA/KT 
CC 26 1 = 1 » K T 
YB=VR+( E! I 1-YBAR 1**2 

26 CCNTINUE 

EMS! X 4- 1 )s=yg/K T 

K T= K T * 1 

URITE<6,3CIK,X(K),Z!K) 9 EMS(KI,U!K) 9 Vm 
30 FORMAT! ?X, 14, 4X, E 13.6, 7X , E 1 3 . 6 , 7X , E 13.6, 7X, E 13.6* 7X, El 3.6 1 
85 CCNTINUE 

PREPARE HEMS FOR PLOTTING 
CC 7 C 1 = 1. N 
T I M fc ( I 1 = 1-1 
7G CCNTINUE 

CC 110 1 = 1, N 
110 ES! I ) = E! I 1**2 
KC=1 
NPT = 300 

DC 43 I=1,NA,10 
YXX! KQ )=EMS( I 1 
43 KC=KC*1 

READ! 5,10111 ITITLE! I). I* 1,121 
101 FORMAT ( 6A 8 1 

CALL DRAWtNPT.TIME.YXX.O.C, 1 
END 



', HI TLE, 0.0, 0.0, 0,0, 0,0, 8, 9, l, LA) 



90 



INITIAL DISTRIBUTION LIST 



No. 



1. Defense Documentation Center 

Cameron Station 
Alexandria, Virginia 22314 

2 . Library- 

Naval Postgraduate School 
Monterey, California 93940 

3. Commandant (PTP) 

U. S. Coast Guard 
Headquarters 

Washington, D. C. 20226 

4. Prof. S. R. Parker (Thesis Advisor) 

Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 93940 

5. Prof. G. S. Thaler 

Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 93940 

6. Prof. H. A. Titus 

Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 93940 

7. LT Harold G. Fletcher, Jr., USCG 

136 Kiwanee Road 

Warwick, Rhode Island 02888 



Copies 

20 

2 

2 

5 

1 

1 

1 



91 



UNCLASSIFIED 

Security Classification 



DOCUMENT CONTROL DATA - R&D 

(Security clameificatton of title, body of abstract and indexing annotation muet 6© entered when the overall report ia c laaeified) 



1. ORIGINATING ACTIVITY (Corporate author) 

Naval Postgraduate School 
Monterey , California 



2a. REPORT SECURITY CLASSIFICATION 

Unclassified 

2b. CROUP 



3- REPORT TITLE 



A Comparative Study of Discrete Time Filtering for a Non-Stationary Random Input 



4- DESCRIPTIVE NOTES (Type ol report and inclueive datea) 

Thesis, Master of Science, December 1967 

5- AUTHORfS; (L set name, li ret name, initial) 

FLETCHER, Harold G., Jr. 



€• REPORT DATE 

December 1967 



la. TOTAL NO. OF PACES 

91 



lb. NO. OF REFS 

19 



8«. CONTRACT OR GRANT NO. 



9a. ORIGINATOR'S REPORT NUMBERfSj 



b. PROJECT NO. 



0 . 



9b. OTHER REPORT NOfSJ (A ny other numbere that may be aeeijned 
thia report) 




11. ABSTRACT 

This thesis is concerned with a comparative study of discrete time filters 
using the theories of Wiener- Kolmogorov, Bode-Shannon, and Kalman, applied to 
the filtering of a non-s tationary random signal in the presence of measurement 
noise. Programs are developed for the simulation of these systems and signals 
on a digital computer. Their filtering properties are compared for a random 
input signal with known steady state characteristics, starting at initial time 
t = t-. The results show that when the gain of the Kalman filter is equal to 
its steady state value, the Kalman and Wiener-Kolmogorov filters perform 
identically. For large initial errors the Kalman filter, with large initial 
gain, gives the best transient response; for small initial errors the Kalman 
and Wiener-Kolmogorov filters are essentially equivalent in their transient 
responses . 



DD 1473 



93 UNCLASSIFIED. 

Security Classification 



UNCL ASSIFIED 

Security Classification 



1 4 

KEY WORDS 


L 1 N 


K A 


LINK b 


LINK C 


ROLE 


W T 


ROLE 


W T 


ROLE 


W T 


13 - 

A * 


4 

4 


- A . 




- *■ 


0 

f ' 

itb 


1 

— 1 
<4 



DD I MOV 65 14 73 (BACK) 94 UNCLASSIFIED 



' t ' I - 0 7 > 6 *» ? 1 



Security Classification 



A 31404 



DUDLEY KNOX LIBRARY 




3 2768 00407336 1 
DUDLEY KNOX LIBRARY 







